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Generation of dark solitons from large initial excitations and their evolution in a quasi-one- 
dimensional Bose-Einstein condensate trapped by a harmonic potential is studied analytically and 
numerically. In the case of a single deep soliton main characteristics of its motion such as a fre- 
quency and amplitude of oscillations are calculated by means of the perturbation theory which in 
the leading order results in a Newtonian dynamics, corrections to which are computed as well. It 
is shown that long-time dynamics of a dark soliton in a generic situation deviates substantially 
from outcomes of the naive application of the Ehrenfest theorem. We also consider three different 
techniques of controllable creation of multi-soliton structures (soliton trains) from large initial ex- 
citations and calculate their initial parameters (depths and velocities) with the use of a generalized 
Bohr-Sommerfeld quantization rule. Multi-soliton effects are discussed. 
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I. INTRODUCTION 



Experimental observations of dark solitons in a cigar- 
shaped Bose-Einstein condensate (BEC) of sodium Q 
and of 87 Rb0 atoms and in a near spherical BEC of 
23 Na atoms Q have stimulated intensive theoretical stud- 
ies devoted to generation and evolution of such excita- 
tions in BEC's. At low enough temperature the conden- 
sate evolution is described sufficiently well by the Gross- 
Pitaevskii (GP) equation p|, which is originally three- 
dimensional (3D) but in some cases is reducible to a ID 
nonlinear Schrodinger (NLS) equation with an external 
potential. As such cases we can mentioned a BEC with 
a "pancake" geometry (see e.g. and a low-density 
BEC in a "cigar-shape" trap (see e.g. @). If the re- 
spective size of the condensate is much bigger than the 
characteristic dimension of an excitation in it, then it is 
commonly believed that the conventional (i.e. without 
potential) NLS equation can be used for understanding 
the excitation evolution during initial intervals of time. 
Such approximation was used for description of dark soli- 
ton dynamics in quasi- ID BEC with a positive scattering 
length (see e.g. |7[). If, however, a size of an excitation 
is comparable with the size of the condensate or if one 
is interested in long-time dynamics, then influence of a 
trap potential must be taken into account. This is es- 
pecially important in the case of dark solitons. One of 
the physical reasons for that is that in the presence of a 
trap potential the BEC density becomes a function of a 
space coordinate and/or time. From the mathematical 
point of view it means that the trap potential changes 
boundary conditions for the macroscopic wave function 
of the condensate at the infinity. (Note that this does 
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not happen in the case of bright solitons in a BEC with 
negative scattering length, which in a quasi- ID case were 
recently observed experimentally |8()- 

In the homogeneous NLS equation, which is exactly 
integrable the long-time evolution of an excitation 
can be reduced to the motion of a set of solitons and 
to propagation of linear waves. Solitons in that case are 
well defined objects and their parameters can be found 
by means of the inverse scattering transform method [jj . 
In particular, constant velocities of solitons (as well as 
other soliton parameters) can be calculated from the as- 
sociated linear spectral problem using the initial distri- 
bution of the macroscopic wave function. In the pres- 
ence of a trap this method cannot be applied without 
some reservations. Moreover, a dark soliton against a 
nonuniform background is not a well defined object any- 
more. Therefore separation of the excitation evolution 
into "soliton part" and "background part" is somewhat 
conditional. Even if such separation is meaningful at 
the initial moment of time, during propagation soliton 
parameters become functions of time. In particular, a 
change of a soliton shape can be so dramatic that the 
distinction between the soliton and the background may 
loose its sense and must be reconsidered. Then one meets 
a problem of long-time description of the soliton evolu- 
tion. This problem becomes even more complicated in 
the case when an initial excitation results in formation 
of many solitons (soliton trains). Thus, one of the aims of 
the present paper is to consider evolution of excitations 
which initially can be classified as one- and multi-soliton 
ones in ID BEC confined by a harmonic potential. 

Several aspects of this problem have already been ad- 
dressed in literature. In particular, there have been re- 
ported several evidences that the dynamics of a dark 
soliton or a few solitons in a BEC trapped in a har- 
monic potential is oscillatory [13, IH El El- 

In one 

of the first publications 10] it has been suggested that 
for description of the dark soliton evolution one can em- 
ploy the Ehrenfest theorem which allows one to define 
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a frequency of oscillations, and on this basis an empiric 
formula m'x — —dU/dx, where U(x) is a trap poten- 
tial, for Newtonian dynamics of a dark soliton has been 
proposed (although a definition for the coordinate of a 
dark soliton x has not been given). A similar equation 
for a dark soliton motion against a background was de- 
rived in [T^ by means of the multi-scale analysis. In 
numerical simulations provided in [13j authors observed 
an oscillatory behavior of small amplitude dark solitons. 
That paper however has left open a question about the 
frequency of the oscillations and their amplitude. Also, 
although an important question about a choice of a par- 
ticular form of the background has been posed (since a 
"wrong" choice results in rapidly growing oscillations of 
the background), it has been solved by numerical means 
on the basis of elimination of the mentioned oscillations 
of the background, leaving analytical fundamentals for 
such a choice open, as well. Moreover in Ref. it has 
been found that the amplitude of a small-depth dark soli- 
ton increases as the soliton approaches a turning point, 
which does not corroborate with the previous findings 
[ToL ITTl IT^ | where it was reported that the shape of a 
dark soliton is preserved during propagation against a 
background. Finally, it is worthwhile to mention that in 
Refs. lid lul lia lK| one can find no information about 
long-time behavior of a dark soliton in a trap potential. 

Another related issue of the theory is generation of 
dark solitons. The discussion in the literature has been 
mainly concerned with such methods as laser induced 
Raman transition between two internal states of the con- 
densate [T3. "phase imprinting" 0, modulational 
instability ||| , and the use collision of two initially sepa- 
rated condensates [Tol.lTT|. 

In Ref. [10( it was supposed that generation of the dark 
solitons with different initial velocities can be achieved 
in a collision of two pieces of initially displaced conden- 
sates. By numerical calculations it was shown that in 
the presence of a trap all solitons have essentially the 
same classical period of oscillations given by harmonic 
classical dynamics. The same mechanism of generation 
of the multi-soliton structures was considered in Ref. [ll| 
(although a number of the generated solitons and their 
initial parameters were not calculated). In Ref. [l4| a 
scheme to create dark solitons in a BEC with the use of 
coherent Raman process which couples internal atomic 
levels of the condensate with a laser was proposed. In a 
recent work on the controlled generation of dark solitons 
by phase imprinting method 15] the authors provided a 
simple theoretical description of the creation of the dark 
solitons in the experiment |2(. Considering a very short 
time scale, they neglected the trap potential, and in this 
case for a given initial phase step, a number of solitons 
and their initial velocities were calculated analytically by 
mapping Zakharov-Shabat eigenvalue problem into the 
pendulum problem, which is mathematically much sim- 
pler. 

Thus, the second aim of the present paper is the study 
of generation of multi-soliton structures from an arbi- 



trary initial excitation in a nonuniform condensate con- 
fined by a harmonic potential as well as their evolution. 

The organization of the paper is as follows. In Section 
II the dynamics of a single dark soliton in a harmonic trap 
potential is considered in detail. First of all we argue the 
choice of the analytical form of the background which 
is necessary for a stable long-time dynamics of a soliton 
making it nearly integrable. We show also that one can 
define two characteristic coordinates associated with the 
dark soliton: a position of the center of mass (the mass 
being negative) and a position of a local minimum of the 
intensity. The both values being the same in the case 
of a homogeneous background display different behavior 
subject to the effect of the potential. The first one is gov- 
erned by the Ehrenfest theorem while another one can be 
described on the basis of the perturbation theory of soli- 
tons ^(| . Theoretical predictions are compared with the 
numerical simulations. Non-adiabatic effects of the soli- 
ton dynamics and the behavior of its phase are discussed. 
In Section III we describe generation of soliton trains in 
a trapped BEC from initial excitations of different types 
and study their evolution. In the case of perturbation 
of the condensate density by a large and smooth ini- 
tial pulse we find initial parameters of created solitons 
with the use of generalized Bohr-Sommerfeld quantiza- 
tion rule |T3. Using the results of Section II we predict 
locations of turning points which are well confirmed by 
direct numerical simulations. We also discuss behavior of 
dark solitons generated by the phase imprinting method 
and creation of solitons during collisions of two conden- 
sates in the presence of the harmonic trap potential. The 
outcomes of the theory are summarized in Conclusion. 
For the sake of convenience a summary of some techni- 
cal results of the perturbation theory for dark solitons is 
given in the Appendix. 

II. MOTION OF A SINGLE DARK SOLITON IN 
A PARABOLIC TRAP 

A. Statement of the problem 

As it is customary, we start with the CP equation for 
the order parameter ip = ip(r, t): 

dih h 2 

ih-^ = -—Ai> + V trap {ryj + g \i;\ 2 i J , (1) 

where we use the standard notations: go = A-Kh 2 a s /m, a s 
being the s-wave scattering length, which is considered 
positive, and m being the atomic mass; Vt rap {r) is a trap 
potential. 

A self-consistent reduction of Eq. (QJ to a ID NLS 
equation can be made in various situations (see e.g. la, 

B). 

(i) In the case of a pancake BEC we suppose that in the 
transverse direction (i.e. in the direction orthogonal to 
the x-axis) the size of the condensate is large enough to be 
considered infinite in the first approximation. This leads 



3 



to the trap potential of the form Vt ra p = ^u^x 2 where 
u)q is the harmonic oscillator frequency. Then, in order to 
rewrite the dynamical equation in a dimensionless form 
we make a substitution 

ip(r, t) = 23 (ATiva s al)~^ exp (ik±r± - ®( x i *)) 

(2) 

where r_L = (y, z) and a 2 , = and make a change of 

independent variables x i— > ^ 1 / 2 aox/2 1 / 4 , i i— > 2 1 / 2 vt/too, 
which results in the canonical form of the NLS equation 
with a parabolic potential 

+ (3) 

(in what follows \& is also referred to as a macroscopic 
wave function). Notice that in contrast to the generally 
accepted renormalization, here we introduced a param- 
eter v which on the one hand characterizes a strength 
of the parabolic potential in dimensionless equation J2J|, 
and on the other hand is connected with the density of 
particles in the condensate. To elucidate this connec- 
tion, we note that the condensate wave function ^(x,t) 
is normalized according to 



/oo 
|*0M)| 2 da; = 4tt2- 1 /V/ 2 
-oo 



a s a n , 



(4) 



where no stands for the "transverse density" of particle, 
i.e. n — N / S, J\f is the total number of particles and S is 
the area of the transverse cross section of the condensate. 
(Formally one has to consider the "thermodynamical" 
limit Af — ► oo, S — * oo at hq =const.) 

On the other hand let us consider an unperturbed con- 
densate wave function ^>(x,t) — exp(— i/it) i i'(x) (/i is a 
renormalized chemical potential) which solves the sta- 
tionary equation 



d 2 * 



2|*| 2 1< = i^x 2 *, 



(5) 



*(0) = 1, *'(0)=0. 



This is clear that decays exponentially at |a;| — > oo, 
and thus one can expect that 



\V{x)\ 2 dx~ / _ \y TF (x)\ 2 dx = i-^ (6) 



3v 



where 



9 TF (x) = \^j2 l i~v 2 x 2 



(7) 



is the condensate wave function in the Thomas-Fermi ap- 
proximation and fj, ~ 2 for More careful numerical 
study of Eq. © gives corrections to © 



\^(x)\ 2 dx 



C(u) 



(8) 



where C(v) is a very slow function of v equal to C(v) ~ 
2.67 with accuracy better than 1% in the interval 0.1 < 
v < 0.3. Comparison of Eqs. J3J and (JHJ) shows that 



2.67- 2 1 / 4 

47T 



2/3 



0.4 



(a s a n ) 2 / 3 (a s a n Q ) 2 / 3 



(9) 



This formula determines v in terms of the physical pa- 
rameters of the system under consideration. As an ex- 
ample one can consider M = 10 5 atoms of 87 Rb (a s as 
5.8 nm) condensed in a trap with ao ~ 1/xm (what corre- 
sponds to the frequency uj ~ 7 ■ 10 2 Hz) and with trans- 
verse radius » 14/im. Then we get v f» 0.2. As it is 
evident, increasing of the density by a factor 10 results 
in decreasing of the small parameter by approximately 
the same factor. 

Thus, v can be viewed as a natural small parameter 
of the problem and therefore in what follows we shall 
consider the case v <C 1. 

In order to construct the theory yet another condition 
is to be imposed: the longitudinal dimension of the con- 
densate ~ v~ x must be sufficiently large compared with 
a characteristic width Z ~ 1 of a typical soliton solution. 
Then it makes sense to speak about a soliton against 
an inhomogeneous background. More specifically, we re- 
quire / /ao ~ y/v. This however is not a strong restriction 
since already at v w 0.1 the relation I « ao/3 (which cor- 
responds to the actual experimental settings) is verified. 

(ii) Considering a cigar-shaped BEC one uses the 
multi-scale expansion (see e.g. 6] for details) with re- 
spect to the small parameter e = 'a s a±/ '<Zq <C 1 (it 
can be viewed as smallness of the energy of two-body 
interaction compared with the kinetic energy) where 
aj_ = (ft/m^jj 1 / 2 is the linear oscillator length in the 
transverse direction and v± is the linear oscillator fre- 
quency in the transverse direction. The smallness of the 
parameter e provides the conditions when the consider- 
ation can be restricted by the lowest mode in the three- 
dimensional parabolic trap: at larger e interactions of 
modes become essential and must be taken into account 
what leads to a set of coupled nonlinear Schrodinger 
equations. In this case the time is measured in units 
2/z/j_ and spatial coordinate along the cigar axis is mea- 
sured in units a±. The parameter v is now defined as 
v = aj_ /ao and can be easily made as small as necessary. 
Meantime, in that case one has to verify the condition 
e<l necessary for applicability of the theory developed 
below. Let us do that for the data analogous to ones re- 
ported in imposing however a lower particle density, 
i.e. considering a BEC of J\f = 5 • 10 5 sodium atoms 
(a s ~ 2.7nm). Then for ao ~ 450 //m and aj_ « 10 /im 
one estimates v « 0.1. Thus in the case at hand the 
limit of small effective parameter v is also reachable for 
available experimental settings. 

In order to complete the statement of the problem we 
have to specify the terminology "long-time" dynamics. 
We will do this for a pancake geometry. A unity of time, 
At = 1, in the dimensionless variables corresponds to 
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ujqIv real seconds. This means, for example, that con- 
sidering a BEC of 87 Rb atoms in a trap with the lon- 
gitudinal size of order of ao ~ lpm and assuming that 
l/ao ~ 0.3, a unity of the dimensionless time will corre- 
spond to 0.18 ms. A typical life-time of a BEC reachable 
in today experiments can be estimated as at least 50 ms, 
which means that it is meaningful (in a BEC context) 
to investigate the dynamics of the excitations for times 
at least up to 300 dimensionless units. As one expects v 
to be of order of the frequency of soliton oscillations in 
a trap, the above estimate corresponds to approximately 
20 periods of oscillations. Thus, long-time dynamics in 
the present paper is understood as a dynamics displaying 
more than 10 oscillations. 



B. Choice of the background and near-integrable 
dynamics 

When v = 0, Eq. reduces to the conventional NLS 
equation which admits a stable constant-background so- 
lution |^(a:,t)| 2 = po(=const), against which one can 
construct a dark soliton solution 9] 



l + ^cxp{ m [x-X(t)]} 
l + exp{rio[x-X(t)]} 



(10) 



Here X(t) = Vt + X a , V = -2^008(19/2) and r] = 
sin(i?/2) are the soliton coordinate, velocity, and 
width, respectively. Xq is the coordinate at initial mo- 
ment of time. Due to the symmetry of the problem, 
the parameter d which characterizes the phase differ- 
ence of fy s (x,t) at ±oo can be restricted to the inter- 
val < < 7r. For d — n the BEC density vanishes 
at the soliton center and in this case a soliton is often 
called "black" while for other $ a soliton is referred to as 
"grey" . The limit corresponds to small amplitude 

solitons. 

When the trap potential is switched on, the function 
$? s (x,t) given by is not a solution of J3J anymore. 
However, one could expect that if a region of soliton lo- 
calization is much less than the size of the background, 
i.e. when v <C 7?0j an d if in the vicinity of the poten- 
tial center the initial condition for Eq. is chosen to 
be close enough to the dark soliton one, then a solution 
of Eq. ||3J) can be searched in a form of a "dark soliton" 
<f>(x,t) against an inhomogeneous background F(x), i.e. 
in the form 



$>{x,t) = F{x) ■ $(x,i), 



(11) 



where $(x, 0) = ^ s (x,0), ^ s (x,t) being given by ifTUI) . 

Ansatz (|llf) is meaningful if the dynamics of &(x, t) is 
close, in some sense, to the dynamics of ^ s (x,t). Since 
^ s (x, t) is governed by the unperturbed (i.e. with v = 0) 
NLS equation, to satisfy the last requirement it is nat- 
ural to choose F(x) such that the resulting equation for 
<S>(x,t) would be close to the NLS equation. This can be 



achieved by requiring F(x) to be an eigenfunction of the 
nonlinear spectral problem [c.f. |J5j] 



F xx + uJb 



{vxf 



F - 2 Po F 3 = 0, 



(12) 



lim F(x) = 0, (13) 

x — *±oo 



which satisfies the following normalization conditions 

F(0) = 1, F x (0) = 0. (14) 

In (|12|l uji, is an eigenvalue. 

Indeed, substituting ansatz ((TTJl in Eq. ®, one obtains 



m + - 2(|$| 2 - p )$ = R[$, F], 



(15) 



where 



R[§, F) = -2 (lnF) x $ x + 2$(|$| 2 - p )(\F\ 2 - 1). (16) 

Next, one can estimate the order of magnitude of R[&, F] 
considering xv small enough: 



(\nF) x ~v 2 x, $ x ~ry, 
(l^-poHlFp-l^W, 



(17) 



where it is taken into account that in a vicinity of the 
center of the potential, i.e. at < 1, F(x) can be 
expanded in the Taylor series 

F(x) = 1 + ^LpV + 0(x 4 ) . (18) 

Substitution of (jl8(l into 112(1 with subsequent expansion 
in powers of "small" x (i.e. x — o(^ -1 )) and small v 
yields the estimate for the eigenvalue 



ui b = 2p + — + 0(is 4 ). 
4po 



(19) 



As it follows from the definition of the soliton width, 
we have ?7 = 0(1) whenever po = 0(1). This means that 
v <C 1 is indeed a small parameter of the problem in 
the sense that R[$,F] = 0(v 2 ). Then the evolution of 
the function ^(x, i) will be described by nearly integrable 
equation 1 (15(1 and in this sense will be sufficiently close 
to the evolution of ^ s (x,t) given by ijTUj) at least for 
period of time defined by t <gC v~ A . We will refer to the 
respective dynamics as nonlinear. If, however, xv > 1 
then the first term in 1(16(1 becomes exponentially small 
(since as it follows from 1(12(1 decay of F(x) at x — > ±00 is 
exponentially fast) and the second term is approximately 
canceled with the nonlinear term in 115(1 . In that case the 
dynamics is governed by the effectively linear equation. 



C. Nonlinear dynamics of a dark soliton. 

Suppose that X(t) — o(^ _1 ). In this case estimates 
((17(1 hold for the whole spatial region and one can sim- 
plify the expressions for the background with the use of 
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and CU): 



F = F (x) + 0(v 4 x 4 ), F (x) = 1 - — x 2 (20) 

8/30 



and rewrite R[&, F] in the form 



F] w F ] = — [* a 



xd>(|$| 2 - Po )] • (21) 



Now Eq. l|15|l can be treated by means of the pertur- 
bation theory for dark solitons jlfij . In this way, however, 
one meets a problem: the term proportional to x$> x in 
the right hand side of l|21[l belongs to the class of pertur- 
bations which effect on the soliton dynamics cannot be 
described by the adiabatic approximation only. To avoid 
this difficulty, we make an ansatz 



i—fft)^- 
2po dx ' 



where 



f(t) = / X(t')dt' + /o, 
Jo 



(22) 



(23) 



and the constant /q is to be determined later. 

Necessity of the renormalization l|22|) can be justified 
by full invoking the perturbation theory (see Ref. 
for details). In order to give here some simple indication 
on the origin of the phenomenon, we notice that secu- 
lar terms appearing due to perturbation of a nonlinear 
equation are eliminated in the so-called adiabatic approx- 
imation which is obtained from the exact soliton solution 
by allowing parameters to be dependent on time. If in 
our case one substitutes e~ l ^/ 2l $> s (x, t), where ^ s (x, t) is 
given by fT7)jl , into the left hand side of 10 and computes 
the imaginary part, one ensures that it is an even function 
of the spatial coordinate. Thus if Im e~ l ®/ 2 R[<&, F] has 
an odd (with respect to x) component, a secular terms 
originated by such a term cannot be eliminated by any 
modification of the soliton parameters. This requires in- 
troducing additional phase factor <p(x, t) [see (|27H below]. 
The peculiarity of such a phase shift, however, is that the 
imaginary part of the respective term (after multiplying 
by exp(— id/2)) is proportional to tanh?7(a; — X{t)) and 
thus is zero at x = X(t). Thus if Ime~ M / 2 R{$, F] is 
not zero at x — X(t), it cannot be eliminated by any 
modification of the adiabatic theory. The ansatz l|22|) 
allows one to count the mentioned "dangerous" term ex- 
plicitly. Taking into account the explicit form of that 
term, namely the fact that it is related to translational 
invariance of the system and is localized about the dark 
soliton kernel, it can be interpreted as an internal mode 
of a dark soliton excited by spatially varying background. 

Now the dynamical equation for <fi with the accuracy 
0(v 4 ) reads 



Po)<t> = v 2 R, 



(24) 



where 
R 



1 



— [{x-X{t))4> x -x\ 



Aif{t)4> 2 4> x 

v 4 f 3 (t) 



v'f 2 (t) 
Po 



(24>\(j) x 



Po)(f> 



2p§ 



4>x\4>x 



(25) 



where the terms proportional to the second and third 
order of fit) in some cases may be not small (see below). 

The perturbation theory can be applied to 
To this end we pass to new independent variables (x,t) — > 
(0,i) where 

e = r)(t)(x-X(t)) and X(t) = Vt + x (t) (26) 

and rj(t) and xo(t) are allowed to be dependent on time 
with the initial conditions 77(0) = 770 and Xo(0) = Xq, 
and look for d> in the form 



(f>ix,t) 



iu <p(x,t) 



V 2 <t>l 



where 



Poe 2 cos 



. . , e 
i sin — tanh — 
2 2 



(27) 



(28) 



is the adiabatic approximation. As it has been mentioned 
above, r\{t) and Xo(t) are (slow) functions of time, de- 
scribing variations of the width/depth and velocity of 
the soliton, respectively. Quantity X(t), which gives a 
position of the minimum of the soliton intensity must be 
interpreted as a coordinate of the soliton center. By intro- 
ducing the rapidly varying phase (pix,t) one can satisfy 
the equation of balance of the momentum |l6j 



P 



1 

2i. 



which reads 



dP 
~dt 



b@<j> - (p<f>e)d@, 



iR(j>e + R<t> @ )dQ 



(29) 



(30) 



(all functions in the integrands are considered to be de- 
pendent on the variables (G,i)). 

In order to find dependence of n(t) on time, one can 
use the conservation law 



where 



dN 
dt 



N = 



dxlmi(/)R) . 



oMpo-H 2 )- 



(31) 



(32) 



Direct calculations allow one to ensure that l|31|l . I|32|) 

give drj/dt — 0, and thus the amplitude of the soliton 
is conserved: r\ = 7/0 • This results gives an explanation 
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of the numerical findings reported in jToL ITTl IT^ | . It is 

important, however, that it has been obtained in the limit 
r] 3> v and thus cannot be compared with the outcomes 
of [3. 

Now, computing R and substituting the result in i|29|) 
and l|30|) . one obtains subject to the boundary conditions 



lim 



|*|- 



e 



4p cosh^(6/2) Jo 



X(t')dt' . 



(33) 



Let us consider a soliton trajectory. After direct sub- 
stitution of Eo.ffiSjl into Eas. (|A2|l . (|A3|I and then substi- 
tution of the obtained result into (|A1|) one can get (see 
Appendix A) 



dxo 
~dt 



— 2 sin^ - + 1 



X{t')dt' 



v 2 V 2 4i/ 4 Fsin 2 f 
— — X 
8po 

2i/ 6 sin 



4 



15p, 







9p 

X{t')dt' 



X(t')dt' 







(34) 



Requiring the initial soliton velocity to be V, i.e. 
^■| t=0 = one obtains the constant introduced in (|23|l : 

- 24 p„^ + <A" K ). 

Below we will argue that the theory is applicable for 
relatively small velocities, in particular for V ~ v. This 
means that we are to restrict the consideration to d close 
to 7r. Let us first consider \d — ir\ ~ v. This allows one 
to verify, using Eq. i|34|) . that X(t) in the leading order 
with respect to v (designated below as X (t)) satisfies 
the equation of the harmonic oscillator 



and thus 



V 

X Q (t) = — sintut) 
v 



(35) 



(36) 



(where the initial conditions Xo(0) = 0, -Xo(O) = V have 
been used). Thus, the soliton center, defined as a local 
minimum of the intensity undergoes oscillatory motion 
with the frequency which is equal to the frequency of 
the trap v. This result also corroborates with previous 
investigations [Hill El- Notice that from Eq. ® and 
the condition |A(£)| < 1 follows the relation V ~ v used 
above. 

Let us now return to expansion (|20H and observe that 
it is still valid if we consider X ~ v^ 1 / 2 (and v 1 ! 2 <C 1). 
In this case the soliton velocity is allowed to be of order 
of j/ 1 / 2 . Then differentiation of l|34|) with respect to time 
yields (instead of 1351) ) 



d 2 X 
~dH r 



v 2 X(t) 



Po 



23 1 

— sin(2z/t) + — sin(3i4) 



(37) 



where 



v 2 



1 



1Y1 
9 ftT 



(38) 



is the renormalized frequency and in the right hand side 
X(t) has been substituted by X (t) given by (J2HJ). Then 
one computes 



X(t) 



1 



541 V 2 \ V 



■ sin(££) 



2160 po 
V 3 ( 23 

— (-sin(2^) + -sin(3^)). (39) 



The obtained result reveals three important features of 
the dynamics when the soliton velocity increases. First, 
the frequency of oscillations increases compared with the 
frequency of the harmonic trap (in the case at hand by 
g^-)- Second, the amplitude of oscillations also increases 
(by the value ^g^)- Third, there exists a small (~ 
v) frequency mismatch between second harmonic and a 
double first harmonic and between third harmonics and 
the triple first harmonic. This leads to slow (compared 
with the period of soliton oscillations) variation of the 
amplitude of the periodic motion. Below we will observe 
all these phenomena in numerical simulations (see Fig.[T]l. 



D. On the Ehrenfest theorem for dark solitons 

As was indicated in jToL IT2]] , the frequency of the soli- 
ton oscillations can be found with help of the Ehrenfest 
theorem which must be written for the "center of gravity" 
of the whole condensate 



x(t) 



1 

N 



where 



N 



x\ty(x,t)\ 2 dx, 



\^{x,t)\ 2 dx 



(40) 



(41) 



is a renormalized number of particles of the whole con- 
densate (N ~ Af) . Then the coordinate of the dark soli- 
ton can be defined as 

1 f°° 

= w x{F(x)} 2 (Mx,t)\ 2 - Po )dx, (42) 



where 



N s = 



(\®(x,t)\ 2 - Po )dx 



(43) 



can be interpreted as a negative mass of a dark soliton. 

One can easily show that x s (t) coincides with Xo(t), 
defined in subsection III CI only in the leading order and 
at initial stages of evolution, and thus cannot be treated 
as a coordinate of the soliton center in a general case. 
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Indeed, it is straightforward to show that x s (t) obeys the 
exact "Newton law" 

f S i =-/I^ |4,(i ' f)|3<fc ' (44) 

which in the case at hand is reduced to the equation for 
the harmonic oscillator. Then, taking into account the 
initial condition in the form 

*{x,0) = F(x)-$ s (x,0), (45) 

which corresponds to the form of the solution given by 
one computes 

V 

x s (t) = — sin(i4), (46) 

which is an exact formula (notice that Eq. I)36|) is approx- 
imate). 

Thus, we have obtained that in the lowest order of v a 
dark soliton with a small velocity V ~ v undergoes oscil- 
latory behavior which emerges both from the perturba- 
tion theory for solitons and from the Ehrenfest theorem. 
Such a behavior is confirmed by the direct numerical sim- 
ulations (see Fig.lc below, where v — 0.2, and V 0.14). 
However, significant deviations from the Ehrenfest theo- 
rem are observed with growth of the initial dark soliton 
velocity [c.f. gBJ and ©]. 



Third, one observes that the amplitude of oscillations 
are larger then V/ v predicted by Ehrenfest theorem and 
displays the periodic modulation. All these differences in- 
crease with growth of the initial soliton velocity V. This 
behavior is also in qualitative agreement with the result 
(|39fl obtained within the framework of the perturbation 
theory. 



FIG. 1: Dark soliton trajectories X(t) (solid lines): for po = 1 
and different initial parameters (a) i?o = 2; (b) i?o = 2.5; (c) 
■do = 3, and the trajectories following from the Ehrenfest 
theorem x s (dashed lines); Figure (d) shows soliton trajecto- 
ries for $o = 2 and different amplitudes po — 0.5; 1; 2 (solid, 
dotted and dashed line correspondingly). The parameter of 
harmonic trap potential is v = 0.2. 

It is to be emphasized that all above effects become 
easily visible after several periods of oscillations, while 
during the first period they closely follow the "Ehrenfest 
trajectory" given by (|46() or l|36[l . This coincidence of the 
trajectories has been observed in Refs. 0, 0, 0, 0] 
where only few cycles were counted. Our numerical sim- 
ulations of first turning point also have shown that if an 
initial soliton is deep enough (what corresponds to small 
values of V), then the amplitude atheor — V/v found from 
Eq. (|36|) and amplitude a num calculated numerically as a 
local minimum of the intensity practically coincide with 
each other, which is illustrated by Table [I] 



E. Numerical simulations of the one-soliton 
dynamics. 

To verify the above findings, we have done extensive 
numerical simulations of the dark soliton dynamics which 
is governed by the evolution equation Eq. subject to 
the initial conditions l|45l) . In Fig. ^ we present a dark 
soliton trajectories X(t), which numerically were defined 
as the coordinate of the local minimum of the soliton, 
from the long-time simulations. Although, strictly speak- 
ing, only the case depicted in Fig. ^ corresponds to 
the perturbation theory developed above, all the plots 
display several features qualitatively coinciding with the 
theoretical predictions. 

First, the choice of the form of the background as an 
eigenfunction of the problem l|12|) indeed allows one to 
follow the dynamics up to hundreds of periods of oscilla- 
tions, if initially soliton had large enough amplitude, i.e. 
one still can identify an entity called dark soliton. 

Second, the dynamics of such defined X(t) displays 
relatively large deviations from that, which would follow 
from the Ehrenfest theorem, when x s is identified with 
the soliton coordinate. Namely, in all the cases the fre- 
quency of the soliton oscillations is higher than the fre- 
quency, predicted from the Ehrenfest theorem (see i|46[) h 
The observed frequency change is of order of a few per- 
cents (~ 3% in Fig. 2J)) which agrees with prediction 
(HP- 



TABLE I: Amplitudes of oscillations calculated theoretically 
atheor — V/u and from the direct simulation of Eq. © 
for po — 1 and different initial values of the parameter #o and 
for two different trap potentials. 



#0 


v — 

Q>theor 


0.2 

Q>num 


V 

Q>theor 


= 0.1 


1 


8.78 


9.17 


17.55 


17.99 


1.25 


8.11 


8.33 


16.22 


16.59 


1.5 


7.32 


7.49 


14.63 


14.96 


1.75 


6.41 


6.56 


12.82 


13.09 


2 


5.40 


5.54 


10.80 


11.01 


2.25 


4.31 


4.42 


8.62 


8.77 


2.5 


3.15 


3.23 


6.31 


6.41 


2.75 


1.95 


1.99 


3.89 


3.95 


3 


0.71 


0.73 


1.42 


1.44 



In Fig. |21 we present snapshots of the dark soliton evo- 
lution for a region of parameters when the perturbation 
theory strictly speaking is not applicable ($o = 2 and 
po = 1 correspond to V 1.08 and thus V > v). One can 
observe rather strong non-adiabatic effects which mani- 
fest themselves in a significant deformation of a soliton 
shape, namely in formation of leading and trailing waves 
before and behind the soliton in the "background" dis- 
tribution. 
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FIG. 2: The reduced shape of the initial distribution com- 
puted as \F(x)\ 2 -\<i>(x,t)\ 2 (solid line) for tf = 2 and v = 0.2 
at (a) t = 0; (b) t = 3; (c) t = 7.3; (d) t = 11; (e) t = 15.2; 
(f) t = 18. Thin dashed lines show the background |F(a;)| 2 . 

Some features of the dark soliton dynamics can be un- 
derstood much better if one investigates in detail the 
turning points. First of all, taking into account that 
only a black soliton has zero velocity, it is natural to 
suppose that at the turning points the soliton becomes 
black. The less is the soliton depth, the farer from the 
center is a turning point, and hence the greater is the 
amplitude of the oscillations. This correlation between 
variations of the depth and of the amplitude becomes ev- 
ident if one compares the results depicted in Fig. [3] and 
in Fig. As it has been mentioned in the previous sec- 
tion, these oscillations are due to mismatch between the 
the first and higher harmonics. In all numerical studies 
carried out we have observed oscillations similar to ones 
shown in Fig.^ They can be explained by the mismatch 
3(V — v). We however did not observe pronounced os- 
cillations due to difference in the second harmonics, i.e. 
due to 2{y — v). We attribute this behavior to the fact 
that non-adiabatic effects due to the third harmonics are 
effectively much stronger than due to the second one. 



FIG. 3: Time dependence of the depth of the dark soliton for 
po = 1, v = 0.2 and different initial parameters $o = 2; 2.5; 3 
(curves 1,2,3 correspondingly). 



FIG. 4: Behavior of the hydrodynamic velocity v = diptot/dx 
in the vicinity of a turning point for initial parameters #o = 2 
and v — 0.3. At is a time counted from the turning point. 
Left column is an output of direct numerical simulations and 
right one is formally obtained form the perturbation theory 
(see the text). 

Finally, an important effect accompanying changes 
of the soliton velocity is a change of the total phase, 
tptot = arg$(x, t), of the solution which is especially sig- 
nificant at the turning points. Our numerical studies of 
the phase evolution are shown in Fig. 0] where in the 
left column we depicted the hydrodynamical velocity v 
defined as v — diptot/dx (see also the next section). As 
one can see, near the turning point hydrodynamic veloc- 
ity changes dramatically making a flip (see Figs. 0Jd,c), 
and after that soliton moves in the opposite direction 
(Fig. 0Ji) . In the right column we show the phase vari- 
ation during the reflection formally computed using for- 



mulas of the perturbation theory, i.e. Eqs. i|10|l . i|22|l . 
(|T7|) . and Comparing the two sets of graphics one 

observes remarkable similarity of the main features (and, 
what is more important, the main scales) of the dynam- 
ics in the both cases. Three differences, however, are 
observable. First, the phase of the direct numerical solu- 
tion displays asymmetry, which is naturally explained by 
the asymmetric non-adiabatic modulation of the soliton 
shape (c.f. Fig.|3J. Second, spatial locations of the turn- 
ing points are slightly different. Third, asymptotic values 
of the hydrodynamical velocity far from the soliton center 
differ as well. These discrepancies are explained by the 
fact that the case under consideration does not satisfy 
the conditions of the applicability of the adiabatic ap- 
proximation, and thus one can expect that the last one 
can give only right orders of magnitude of the parame- 
ters and indicate correctly the main qualitative features 
(what is indeed observed in comparison of the two sets 
of graphics). 

III. FORMATION OF SOLITONS FROM LARGE 
INITIAL EXCITATION 

The problem of evaluation of parameters of dark soli- 
tons formed from a large initial excitation on a constant 
background is formally solved by the inverse scatter- 
ing method 9]. In the framework of this method, the 
NLS equation is associated with the so-called Zakharov- 
Shabat linear spectral problem and soliton param- 
eters are related with the eigenvalues of this problem 
calculated for a given initial condensate wave function 
^(x, 0). If the initial disturbance is large enough, so that 
the linear spectral problem possesses many eigenvalues, 
then a well-known quasi-classical method can be applied 
for their calculation. As was shown in recent paper [l7| . 
a generalized Bohr-Sommerfeld quantization rule is very 
convenient for this aim. 

To formulate this rule, it is convenient to introduce a 
new small parameter e, e -C 1, into Eq. by means 
of replacements x — x' /e, t = t' /s, v — v's, so that the 
equation transforms to 

-1 +41-21*!'* = ^*, (47) 

where we have omitted for simplicity the primes in the 
new variables. Then the limit e <C 1 corresponds to 
formation of a large number of solitons from an initial 
disturbance with parameters of order of magnitude 0(1) 
(see |17|). In framework of the inverse scattering trans- 
form method the NLS equation, i.e. Eq. I|47|) with the 
zero right hand side, is treated as a compatibility con- 
dition of two linear equations for auxiliary function \, 
which we write down in the form 

e 2 Xxx = Ax, Xt = -TjBx* + &xX, (48) 
(equivalent to the Zakharov-Shabat problem; see 0, 
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IT^l. where 



A 



B = 2A 



2* 



+ 1*1 



2* 



(49) 
(50) 



The first equation l|48|l may be considered as a second 
order scalar spectral problem with a given "potential" ^ 
and A playing the role of the spectral parameter. When 
*f>(x,t) evolves according to Eq. I|47|l with v = 0, the 
eigenvalues A„ of this problem do not change with time 
t, and each eigenvalue corresponds to a soliton created 
from the initial pulse. To investigate the limit e -C 1, let 
us represent the condensate wave function in the form 



V(x,t) = \/pOM)exp ( -_ I v(x',t)dx'), (51) 



where p(x, t) has a meaning of the condensate density and 
v(x, t) is the hydrodynamic velocity. Indeed, substitution 
of ifST]) into (|4*7Jl yields the system 



\pt + {pv)x = 0, 

\v t + vv x + p x +e 2 ((pi 



2pp xx )/8p 2 



\v 2 x, 



(52) 

which for e — > takes the form of hydrodynamic equa- 
tions. For smooth enough functions p{x,t) and v(x,t), 
when 



\ep x /p\ < p, and \ev x \ < p, 



(53) 



what corresponds to neglecting the space derivatives in 
A, spectral problem (|48|l transforms into 



£ \xx — 



- (A + v/2) 2 + p 



(54) 



It is to be mentioned here that conditions (|53[) are noth- 
ing but the conditions of applicability of the well-known 
hydrodynamical approach when due to a relatively high 
density the two-body interactions are strong enough and 
one can neglect the "quantum pressure" (see, e.g. 0). 



FIG. 5: Schematic plot of the Riemann invariant A + given by 
Eg. 1551 (thick solid line); thin solid line shows the Riemann 
invariant for disturbance with respect to uniform background; 
dashed line shows background without disturbance, x* de- 
notes the position of localized disturbance and AF(x*) — 
1 — F(x*) is the change of the condensate density because 
of condensate non-uniformity, d is characteristic scale of the 
initial disturbance. The horizontal lines of different width in- 
dicate positions of eigenvalues \ n and the width of each line 
characterizes the lifetime of the corresponding solitons (the 
thicker is a line, the smaller is the lifetime). 

Equation (|54|l has a formal analogy with a stationary 
Schrodinger equation for a quantum-mechanical motion 
of a particle in the "energy-dependent" potential, i.e. the 



potential depending on the spectral parameter A. Ac- 
cording the mentioned above independence of the eigen- 
values A„ on time, they can be calculated with the use of 
the initial distributions p(x,0) and v(x, 0). Since we are 
interested in such initial data which give rise to creation 
of large number of solitons, this means that functions 
p(x,0) and v(x,0) correspond to the problem (|54|l with 
a large number of eigenvalues. Then the quas i-classical 
approach can be used for their calculation [13. To clar- 
ify this method, we have shown schematically in Fig. [31 a 
plot of the "Riemann invariant" 



A+ = -v(x,0)/2+y/p(x,0) 



(55) 



which plays a role of the "quantum-mechanical poten- 
tial" for the problem 1)54(1. (The second Riemann invari- 
ant A" = -v(x, 0)/2 - y/p{x,0) can be considered in the 
same way.) The Riemann invariant for the same distur- 
bance but with respect to uniform background F(x) — 1 
is shown for comparison by thin solid line. We see 
that in both cases there is a "potential well" , but for 
the nonuniform background case the eigenvalues acquire 
imaginary part ("decay width") due to tunnelling effect. 
This means that dark solitons in confined condensate 
have a finite life-time r„ determined by the imaginary 
part of the eigenvalue A„ , what correlates with the above 
mentioned fact that in this case the soliton is not a "well- 
defined" object. Nevertheless, it makes sense to speak 
about solitons in a confined condensate, if their life-times 
T n are much greater than the period ~ 2tt jv of their os- 
cillations discussed in the preceding Section. It is clear 
that "shallow" solitons with small r„ do not survive in the 
confined condensate, so that A n with values close to the 
top of the "potential barrier" do not correspond to any 
real solitons. On the contrary, in the case of the uniform 
background discussed in |17| all eigenvalues correspond to 
real solitons appearing eventually from the initial pulse. 
In multi-soliton problem there is also one more scale of 
time equal to time of formation of solitons from the ini- 
tial pulse. For deep enough solitons it can be estimated 
by the order of magnitude as time necessary for solitons 
with velocity \V n \ — 2|A„| (see below) to pass the dis- 
tance equal to the width d of the initial problem. For 
the problems under consideration, this time ~ rf/2|A n | 
must be much less than the period 2it/i>. Thus, we are 
interested in the eigenvalues A n which satisfy inequalities 



2tt d 
r n » — » — - 
v 2A„ 



(56) 



It is clear that there are no these reservations in the case 
of uniform background [T^ | where formally r n — 00 and 
one can wait long enough to observe formation of soliton 
with any value of A n . 

To calculate the real parts of the eigenvalues A„ cor- 
responding to deep solitons, one can use the generalized 
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Bohr- Sommerf eld quantization rule 




p(x, 0) dx = 7T I n + 



» = 0,1,2,. 



,M, 



(57) 



with given initial distributions /5(x,0) and v(x,0). We 
suppose here that the integrand has only one maximum 
and x^ and x~ are the points where the integrand func- 
tion vanishes: they depend on A n and are chosen so that 
the relationship l|57l) is satisfied (see Fig.[SJ). Analytical 
form of each emerging soliton in an asymptotic region 
where it is well separated from other solitons (i.e. in the 
limit t — ► oo) is expressed in terms of A„ as follows: 



P^(x,t)=p 



Po 



A 2 



cosh yjp - \ n (x - 2X n t) je 



> { r ) {x,t) = \ n ( K p /p^ 



- 1 



> (58) 



(59) 



As it is clear, formulas f5%|) . f5§|) and ifTHjl with ■d = d n , 
f]0 = In, and V = V n , where the parameters are con- 
nected by the relation A„ = ^/po cos($ n /2), represent 
the same one-soliton solution. Thus, the last formula al- 
lows one to find initial values of d n for solitons emerging 
from the dark excitation of the condensate with given ini- 
tial distributions of p(x, 0) and v(x, 0) against a constant 
uniform background. 

If the background is not uniform but changes in space 
in the intervals of integration in (|57|l , then we can apply 
the same method with po replaced by the value of the 
background density F 2 (x*) at the place x* of the local- 
ized initial excitation (see Fig. 0. 

We have used this approach for finding soliton param- 
eters for different types of initial excitations: (i) exci- 
tations of the density p(x) (ii) excitation of the hy- 
drodynamic velocity v(x) ("phase imprinting method") 
0, E| and (iii) collision of condensates • 

(i) In the first case of the density disturbance the initial 
data were taken in the form 



p(x,0) 



1 



cosh(x) 



v(x,0) = 0, 



(60) 



where the parameter a measures the strength of the dis- 
turbance. We have chosen the following values of the 
parameters: a = 0.8, v = 0.2, e = 0.3. The values 
of A n for the three deepest solitons calculated with the 
use the Bohr-Sommerfeld rule l|57|) are shown in Table ITT1 
together with the corresponding values of $ n and am- 
plitudes aj^Lr calculated for these soliton from (|36*|) and 
found from numerical solution of @ with the initial data 
iJBSJ. The discrepancy is less than 10% and is caused, ap- 
parently, by the fact that in our case the created solitons 
did not reach yet the asymptotic values of their veloci- 
ties V n = 2A„. Besides that, numerical calculations show 



TABLE II: Parameters of solitons created from initial inten- 



sity disturbance 


n 


A (») 




(n) 
a theor 


(n) 





0.41 


2.29 


2.74 


2.48 


1 


0.65 


1.72 


4.35 


4.23 


2 


0.80 


1.28 


5.42 


5.54 



that the "initial coordinates" of solitons created from the 
initial pulse cannot be identified exactly with X(0) = 0. 
This is the reason why solitons created from one initial 
pulse do not reverse simultaneously their directions of 
motion even during the first period of oscillations in the 
confined condensate. This phenomenon is illustrated in 
Fig. H3 where the density p(x, t) and hydro dynamic ve- 
locity v(x,t) of the condensate are shown as functions of 
x at the moments when two solitons in each train have 
already reversed the direction of their propagation and 
are moving to the center, while the other two solitons are 
still moving from the center. This process of "solitons re- 
flection from the potential well" takes about 20% of the 
whole period 2n/v of their oscillations. 



FIG. 6: Space distribution of the density of a BEC (a) and of 
its hydrodynamic velocity (b) in the harmonic trapped poten- 
tial with v = 0.3 at time t = 4.24 with initial excitation taken 
in the form of pulse 16011 with po — 1, a = 0.8 and e = 0.3. 



FIG. 7: Space distribution of density of BEC (a) and its hy- 
drodynamic velocity (b) in the harmonic trapped potential 
with v = 0.3 at time t = 1 without (dotted line) and with 
(solid line) initial excitation that is taken as a phase step 
with po = 1, a = 4, k = 0.4 and e = 1. 

(ii) For illustration of the process of solitons forma- 
tion by the phase imprinting method we have chosen the 
initial conditions in the form 



p(x,0) = F 2 (x), v(x,0) 



cosh (kx) 



(61) 



Again the solitons parameters can be calculated with the 
use of the Bohr-Sommerfeld quantization rule (|57|> and 
their values correspond well to numerical simulation. In 
particular, the Bohr-Sommerfeld rule gives correct num- 
ber of solitons and signs of their initial velocities. As 
one can see in Fig. [7| at first the deepest soliton moves 
to the right but the hydrodynamic velocity distribution 
corresponds to its motion to the left. Only when the 
local density minimum touches the x axis, the velocity 
distribution makes a flip and after that it corresponds 
to the predicted direction of the soliton propagation (see 
Fig. EJ. All solitons predicted by the Bohr-Sommerfeld 
quantization rule can be observed during some time in- 
terval after their formation (see Figs. I8l9|l . However, in 
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FIG. 8: The same as on Fig. □ at t = 1.5. 



FIG. 9: The same as on Fig. □ at t = 4. 



the case of initial data (|61|) we observe also strong non- 
soliton contribution into excitation of the condensate (see 
Fig. EJ which leads to much more complicated picture of 
its evolution. One may say that in this case solitons move 
along background varying with time and the influence of 
this time dependence is not small on the contrary to the 
previous case of the density disturbance. As a result, 
the motion of solitons cannot be described as (almost) 
harmonic oscillations along constant nonuniform back- 
ground. But even in this case quasi-classical method of 
finding the solitons parameters yields quite accurate de- 
scription of solitons motion for times less than 2-k/v. 

(iii) Formation of dark solitons can be observed also 
during collision of two separated condensates which move 
under influence of the trap potential [Tl|. Taking the 
initial conditions in the form 

p(x,0) = exp \—k(x - i) 2 \ + cxp [—k(x + £) 2 ] , 
v(x,0) = (62) 

we have solved Eq. @ numerically. It was found that 
these two condensates oscillate in the trap potential and 
interact with each other in quite complicated way when 
they overlap with each other. Since the initial data Q62|l 
lead to a number of eigenvalues A*-™- 1 , one may expect that 
during the collision of two condensates the corresponding 
number of dark solitons must be observed. This is indeed 
the case as one can see in Fig. where the density 
and the hydrodynamic velocity distributions are shown 
at the moment of maximal overlap of two condensates 
whose initial density distributions are indicated by dotted 
lines. The number of solitons matches very well with that 
predicted by the Bohr-Sommerfeld quantization rule, but 
their motion cannot be presented as propagation with 
slowly changing parameters along constant nonuniform 
background. 



FIG. 10: Formation of dark solitons (solid line) during the 
collision of two initially separated condensates (dotted line) 
corresponding to the initial data 1621 with parameters n — 
0.5, £ = 8 and e = 0.5 in harmonic trap potential v — 0.2 at 
t — 4. a) density of BEC; b) hydrodynamic velocity. 

Thus, quasi-classical approach provides simple and ef- 
fective method of calculation of the parameters of soli- 
tons arising from large enough initial disturbance. This 
method can be used for estimation of these parameters 
in the today experiments with BEC solitons. 



IV. DISCUSSION AND CONCLUSION 

In the present paper we have investigated the evolution 
of localized excitations in a quasi- ID BEC with a positive 
scattering length confined by a harmonic trap potential. 
It has been shown that in the case of a single dark soliton 
the evolution can be considered as a newtonian one only 
at very low velocities and, hence, big depths of the soli- 
ton and at a large enough effective longitudinal size of the 
condensate. Then the dynamics becomes near-integrable 
and the perturbation theory for dark solitons can be 
used. Non-adiabatic effects become essential when one 
considers long-time dynamics, i.e. dynamics when soli- 
ton makes several oscillations. These effects are: increase 
of the frequency of soliton oscillations; increase of the am- 
plitude of oscillations, which, besides that, changes peri- 
odically with a frequency much less than the frequency 
of the oscillations themselves; change of a soliton shape 
during the evolution. 

Main effects observed in non-newtonian dynamics of 
a soliton are described qualitatively by the perturbation 
theory for dark solitons. In particular, the theory allows 
one to justify the choice of the shape of the background 
which supports many-cycle dynamics of a soliton with- 
out substantial change of its shape. The present study, 
however, leaves open a question about non-adiabatic de- 
formation of a soliton shape. This requires careful study 
of the first order perturbation theory which can be made, 
say, with the use the Green function approach (see 0) 
for the respective linearized problem. 

The existence of an inhomogeneous background be- 
comes especially important when initially multi-soliton 
pulses are under consideration. Then, compared to the 
integrable case with constant background, new tempo- 
ral scales appear in the problem. They are associated 
with the harmonic oscillator frequency and finite life- 
time of solitons. The life-time decreases with the soli- 
ton amplitude, which leads to rather rapid disappearance 
of shallow dark solitons. Generally speaking, a soliton 
with a small amplitude can even loose its meaning at all 
when it is considered against a nonuniform background. 
Indeed, returning to the condition of smallness of the 
soliton width ~ l/rj compared with the linear oscillator 
length ~ l/i/ , one must require 77 3> v for a soliton to be 
meaningful. Recalling now that a t heor < ^/po/v where 
o-theor is an amplitude of soliton oscillations, what fol- 
lows from the Ehrcnfcst theorem, one concludes that it 
make sense to speak about the oscillations of a soliton in 
the case when atheor < 1/^, i- e. when po ~ 1. If po is 
large enough, then the small amplitude soliton can reach 
a region of an exponentially decaying tail of the conden- 
sate distribution already during the first half-a-period of 
oscillations. In that region the dynamics is essentially 
linear and thus the pulse will disappear due to disper- 
sion effects as it happens with linear wave-packets. In 
order to estimate amplitudes of the background at which 
such behavior is observable, we take into account that the 
effective nonlinearity, determined from i|15|) . can be esti- 
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2„2 
die 



and thus it become of order 



mated as po exp(— v 2 d: 
of 10~ 2 already at p ~ 6.5. 

In the context of the above findings a natural ques- 
tion arises about detecting soliton parameters. In con- 
nection with this question it is relevant to mention that 
motion of the soliton is accompanied with the hydrody- 
namic flow with velocity v(x, t) and the corresponding 
matter current, which density in the dimensionless units 
is J(x, t) — p(x, t)v(x, t). Dependence of these two quan- 
tities on time is periodic. Their dependence on the spatial 
coordinate is also non-monotonic. Thus, by detection of 
the velocity and current distribution allows one to make 
an estimate of the dark soliton parameters. 
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APPENDIX A 



For the sake of completeness here we reproduce the dy- 
namical equation for xq (t) which is derived in 0] (notice 
that the sign between the second term in the right hand 
side, is corrected): 



dx 
dt 



dQl R'+(0,t) 

I VoV 

( ©/2 



-5- I — ~o"n + tanh ^ ) R" (6, t) 
V 2 oV I cosh 2 (f ) 1 1 



Vo Jo 



V 2 



sinh 8 
3 1 



6 1- 



4p V' 2 cosh 2 (f) 




, v 2 \ . e 

1 tanh — 



Po 



where notations 



R'±(e,t) = -Re\eR(Q,t)±eR(-G,t) 



R'l(e,t) 



— Im 

2 



eR(Q,t)TsR(-@,t) 



(Al) 

(A2) 
(A3) 



Here R is determined by Ea. (|25ll . 
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